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Abstract 

We apply the Broad Histogram Method to an Ising system in the context of the 
recently reformulated Generalized Thermostatistics, and we claim it to be a very 
efficient simulation tool for this non-extensive statistics. Results are obtained for the 
nearest-neighbour version of the Ising model for a range of values of the q parameter 
of Generalized Thermostatistics. We found an evidence that the 2D-Ising model does 
not undergo phase transitions at finite temperatures except for the extensive case 
Q = l. 
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1 Introduction 



Many systems seem to be well described by a non-extensive thermostatistic 
rather than the usual Boltzmann-Gibbs statistics (BGS). When the effective 
microscopic interactions and microscopic memory are short-ranged and the 
boundary conditions are non-(multi)fractal [1], the BGS provides a complete 
and consistent description of the system. Otherwise, it fails. An alternative 
approach is the use of the so-called Tsallis statistics (TS) [2,3]. The entropy 
in the TS is defined as [2] 



q-l 



with X^iPi = 1) where ? is a given state with energy e, from Vt possible states 
and A; is a positive constant. The index q characterizes the degree of nonexten- 
sivity. The limiting case g = 1 recovers the usual BGS entropy definition. The 
other contraints needed to obtain the thermodynamical averages have been 
recently discussed by Tsallis, Mendes and Plastino [4]. 

Magnetic systems [5]-[19] are to be counted among the large variety of systems 
to which TS has already been applied [3]. The reason for such an interest is 
obvious: statistical physicists have spent a lot of time in the past developing 
tools to understand the critical phenomena presented by this kind of systems 
- and have succeeded in this task. As a natural consequence, one could expect 
some attempts at generalizing well established approaches developed for ex- 
tensive systems. The real space renormalization group [9] and the mean field 
approximation [15] are examples of tools that have been considered in these 
generalizations. However, it is not surprising that results from these approxi- 
mations seem to be controversial: even fundamental questions such as how to 
define and to obtain the correct expression for the thermodynamical averages 
have been revisited very recently[4]. 

The Monte Carlo method is a powerful tool frequently used in the BGS frame- 
work in order to solve ambiguities of this sort. However, within the TS frame- 
work, Monte Carlo simulations of magnetic systems have not been exploited 
so far because of technical difficulties that we are going to address in this 
work. One issue is how to the define the acceptance probabilities which lead 
to the correct distribution of probabilities. A first alternative was presented in 
the application of the generalized simulated annealing for the Traveling Sales- 
man Problem[21], where the acceptance probability is a simple generalization 
of the Metropolis algorithm for q ^ 1. This alternative was motivated by the 
definition of the thermodynamical averages in use at that time, which imposed 
an undesirable dependence on the definition of the zero of the energy scale. 
The use of the above mentioned alternative for the acceptance probabilities 



conveniently removes this dependence. A second alternative was proposed by 
Andricioaei and Straub (AS) [19] a couple of years after the first one. A new 
expression for the acceptance probabilities was obtained from the detailed bal- 
ance condition (a sufficient but not necessary condition for thermodynamical 
equilibrium). This new alternative cleverly circumvents the ambiguity on the 
definition of the lowest level of energy. 

The second reason why, in our opinion, Monte Carlo simulations are not fre- 
quently used within the TS is the fact that all the computational effort in- 
volved in a BGS simulation has to be spent again for each value of the param- 
eter q. Since g is a continuous variable, computer simulations are much more 
time-consuming in the TS than in the BGS, if one wants to investigate the q 
dependence of the results. That is why it has been difficult, for instance, to 
answer the question about which acceptance probabilities to use in a Monte 
Carlo simulation of magnetic systems. One of the goals of the present work is 
to show that the recently proposed Broad Histogram Method (BHM) [22] is 
the ideal tool for Monte Carlo simulations within the TS framework. Because 
the g-independent density of states g{E) is, in the BHM, directly obtained 
from some measured microscopic quantities, we are able to obtain any ther- 
modynamic observable for all values of q and T from only one computer run. 
This fact turns BHM simulations on the TS , for all values of q, as fast as in 
the BGS. 

Briefly, in this paper we are suggesting the BHM as the ideal tool for computer 
studies on the TS. We show it through a simulation of the two-dimensional 
Ising model with short-range interactions. We are aware of the fact that the 
model we are going to study is an extensive one, therefore well described by 
the BGS. Among the reasons why we decided to study this system are: 

• this model can be very easily simulated with great efficiency; 

• the exact solution for the density of states of finite systems is known in the 
limit q = 1 [25]; moreover, the BHM is able to reproduce this exact solution 
with great accuracy [22]; 

• previous results for this model using other approximation methods are con- 
troversial and/or inconclusive; our simulations could shed some light on this 
ongoing discussion; 

• we could easily and reliably show which choice for the acceptance probabil- 
ities reproduces the Tsallis distribution of statistical weights [21,19]; 

• this simple system is here used as a testing ground for the methods we 
propose; building on this first step allows us to use the same approach to 
study , other more complex non-extensive systems such as the Ising model 
with long-range interactions. 

In summary, our choice of a well known system is most convenient since we 
are dealing with two brand-new recipes: the TS with normalized q-expectation 



values and the BHM. 

This paper is organized as follows: in the next section, we review the TS with 
normalized q-expectation values. We also include a discussion about the sta- 
bility of the solutions for the free energy in this new formalism. In section 3, 
we review the Broad Histogram Monte Carlo Method and present its imple- 
mentation for the TS. This is followed by a presentation of our results and 
conclusions. 



2 The Tsallis Statistics with "normalized q-expectation values" 



Tsallis, Mendes and Plastino [4] have recently discussed the role of constraints 
within TS. In that work, they study three different alternatives for the internal 
energy constraint. The first two choices correspond to the ones which have been 
applied to many different systems in the last years [3]. They are: (ia) Yl,iPi^i = 
U and (ib) J^iPf^i = Uq. However, both constraints present difficulties. A third 
choice for the internal energy constraint is defined as [4] 



U. 



l^i=l Pi 



(2) 



where q is the degree of non-extensivity, also present in the definition of en- 
tropy (eq. (1)). Each constraint (ia),(ib) and eq.(2) determines a different set 
of probabilities pi for each state with energy Cj. The extremization of the gen- 
eralized entropy (1), under constraint (2) gives us an implicit equation for the 
probabilities pf. 
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(4) 



The normalized q-expectation value of an observable is therefore defined as 

l^i=iPi^i 



o, = m. 



^ ^1 
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(5) 



where O is any observable which commutes with the Hamiltonian - otherwise 
we should make use of the density operator p. We will refer to this reformula- 
tion of the TS as "with normalized q-expectation values" . A very important 
consequence of this new definition of constraints is that the probabilities do 
not depend on the choice of the zero of energy. 

In order to solve eq. (3) Tsallis et al. suggest two different approaches, namely 
the Iterative Procedure and the (3 ^ 13' transformation. In the iterative pro- 
cedure, we start with an initial set of probabilities and iterate them self- 
consistently until the desirable precision is reached. In the (3 ^ (3' transfor- 
mation the set of equations above is transformed to: 

p,= [l-(l-g)/3'e,]^/Z; (6) 

^; = E[l-(l-g)/?'e.]^ (7) 



with 
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This set of equations is similar to the one that is obtained using constraint 
(ib), except for its dependency on the renormalized temperature, given by eq. 
(8). 

In order to obtain pj, we go through the following steps: 

(1) Compute the quantities |/j = (1 — (1 — g)/?'ej), \/i G Vt] 

(2) If i/i < them Vi = 0; 

(3) Compute Z,,=j:?=iyl^^'~'^ ■ 



g'' 



(4) Compute Pi{p') = y]'''^"'^ /Z, 

(5) Obtain Uq{l3') and any other thermodynamical quantity using eq. (5); 

(6) Obtain /3(/3') from equation (8). 

This recipe allows the determination oi pi{j3) for all /3(/5') and consequently 
Ug{f3) (and any other observable). The second step in the above procedure is 
the well known cutoff [4] associated to "vanishing probabilities". This cutoff 
is required only for g < 1. Because the cutoff is applied before the actual 
computation of the probabilities, the norm constraint is still respected. 

It has been shown [26] that both recipes, the iterative method and (3 —^ f3' 
transformation, demand a careful analysis before its application. The free en- 
ergy obtained from the iterative procedure presents a non-physical disconti- 
nuity whereas the free energy from the j3 ^ (3' transformation has loops. To 



get rid of these pathologies, one has to make exphcit use of the miniinization 
condition on the free energy, whenever an ambiguity appears. This procedure 
generates the correct internal energy and temperature dependency and re- 
stores the proper behaviour of the thermodynamic observables. Following the 
suggestion of ref . [26] , we use in this paper the P ^ P' transformation with the 
proper corrections. 



3 The Broad Histogram Monte Carlo Method 



The approach that we are going to discuss in this section, the Broad Histogram 
Method (BHM) [22]- [24], is one of the many attempts at doing very efficient 
simulations. In traditional simulations, we need a new run for each value of the 
temperature. However, some different approaches have been proposed in which 
we compute some quantities at a given temperature and reweight them for a 
different temperature (see, for instance, ref. [27] and references therein). One 
of these approaches is the histogram method, first introduced by Salzburg [28] 
and popularized by Ferrenberg and Swendsen [29]. However, it has been shown 
[30] that the histogram method presents some limitations, the most important 
concerning the range of temperatures for which just one run is sufficient. The 
BHM enables us to directly calculate the energy spectrum g{E), without any 
need for a particular choice of the thermostatistics to be used [22]- [24]. 

In the BHM the energy degeneracy is calculated through the steps: 

Step 1: Choice of a micro reversible protocol of allowed movements in the state 
space of the system such that changing from an Xom to an X^^w configuration 
is allowed if, and only if, the reverse change is also allowed (the protocol must 
be micro-reversible): 

^ ^ ^ ^ ^ ^ 

allowed allowed 



it is important to note that these movements are virtual, since they are not 
actually performed. 

Step 2: Choice of a fixed amount of energy change AEf^^ and computation of 
A''up(X) (A^dn(^)) for the configuration X, defined as the number of allowed 
movements that would increase (decrease) the energy of the configuration by 
Ai?fix. Then {N^p{E)) ((A^dn(-E))) is the micro canonical average of A^up(-^) 
(N^niX)) at energy E; 

Step 3: Since the total number of possible movements from level E + AE^yr to 
level E is equal to the total number of possible movements from level E to 



level E + AEfix, we can write down the equation 

g{E){N^^{E)) = g{E + AEii,){Nd,,{E + AEf,,)). (10) 



This relation is exact for any statistical model or any energy spectrum [24]. It 
can be rewritten as 

In g( E + AEii,) - In g(E) = In- (^"p(-^)) (n) 

^^ ^ ^^ ' {NdniE + AEf,,)) ^ ^ 



If we choose Ai^fix ^ E, the above equation can be approximated by 

dlng(E) _ 1 ^^ (iVup(i^)) .^2) 

dE AEfix {Ndn{E)) ^ ' 



This equation can be easily solved for g(E). Once this quantity is known, the 
expected value of some observable O can be calculated by 

_ EeJPMemeji 



This method (in the q=l BGS) was first applied to systems with discrete 
degrees of freedom. Recently it has been extended to continuous systems, such 
as the XY Model [31]. To our knowledge, this is the first time the method is 
being applied to a different statistics. Besides the more accurate and faster 
results in comparison to traditional methods, the BHM is even more efficient 
in this particular case because eq. (13) is the only quantity to be recalculated 
for each new value of q. 



4 Implementation of the BHM on the 2D Ising Model with first 
neighbour interactions 



To clarify the application of our ideas to magnetic systems, we will use the 
square lattice ferromagnetic Ising Model with first neighbour interaction. The 
Hamiltonian for this system is given by 

n = H/J=-J2aiaj (14) 

where ctj = ±1, J is a positive constant. The sum is performed over all pairs 
of first neighbours in a square lattice of size N = L x L. For an efficient 



implementation, we rewrite the Hamiltonian as 

Hm = HJJ = J2Q^Q= ^'^ : ' (15) 

{ij) ^ 



where Ci = or 1 and ® represents the exclusive OR operation, where ® 

= 1 (S> 1 = and 1®0 = 0(S>1 = 1. The critical temperature for q = 

1 in the thermodynamical limit for this renormalized Hamiltonian is Tc = 
[2arctanh( V2 - 1)]"^ = 1.13459.... 

We choose a single spin flip protocol of movements to obtain {N^p{E)) and 
{N^^{E)). As proposed in [22] and [23], an unbiased random walk is them per- 
formed on the energy axis in order to visit all values on the energy spectrum. 

The number N^p{E) {Ndn{E)) of movements that increase (decrease) the cur- 
rent value of energy E by AEfi^ is used to calculate g{E). The magnetization 
M{E) is also stored for each lattice size. Here we choose AEf^^ = 4 out of 
the possible values |Ai5|=0,2,4. The histogram is obtained from 6.400.000 
samples (distributed in the energy axis) for L=30,50,70,100. Obtaining the 
histograms for L = 30 takes 20 minutes of CPU time on a DEC Alpha 400. 
For L = 100, the time increases to 220 minutes. 

The next step is to obtain the set of probabihties using the P ^ P' procedure. 
After this step, we can use eq. (13) to obtain the internal energy Uq{T) and 
the magnetization Mq{T). The free energy Fq{T) is also calculated using 

1 Z^-'i - 1 

Fq ^ Uq-TSq = Uq - ^ ^^^ (16) 

Instead of using eq. (4), it is more efficient to use the relation [4] 



1=1 E 

The CPU time for the determination of the values of any observable in the 
whole range of temperatures is independent of the lattice size. Typically, it 
takes 30s on a DEC Alpha 400 for 20000 values of temperature, for each value 
of q. 

For the sake of comparison, we have implemented multispin [32] versions of 
the Ising Model with the AS acceptance probability [19] 

p = ^[1 - tanh{p'AE/2)] (17) 



where E = [1/P'{q 
probability [21] 



1)] ln[l — (1 — q)P'E], and with the TSP acceptance 



P 



[l-q)/3'AE)- 



Notice that the equations above are written as functions of j3' instead of (3. 
The reason is that they were proposed before the pubhcation of the solution 
for the constraints problem in the TS [4] discussed in section 2 . 

For an L = 30 lattice, a simulation run, using the traditional method with the 
same numer of Monte Carlo steps as used in the BHM simulation, took 60s of 
CPU time for a single pair of values (g, T) . This should be compared to 
the 1200s of BHM for all values of q and the whole range of temperatures. 



5 Results 



We will now discuss the results obtained from the implementation described 
in the previous section. We begin with a comparison between the results ob- 
tained through the use of the two already mentioned techniques: the BHM 
and the traditional multispin simulation with both the AS and the TSP prob- 
abilities. The comparison must be made in terms of T', because the AS and 
TSP probabilities were derived before the publication of Ref. [4], as pointed 
above. In Fig. 1 we show the results for the magnetization and the internal 
energy as functions of T'. 
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Fig. 1. Comparison between the results from the BHM, AS and TSP approaches. In 
both the AS and TSP simulations, 1000 samples were used for each T. The BHM 
and AS probabilities agree in both q < 1 and q > 1 regime. The TSP probabilities 
deviate from these results even for values of q very close to one (where the three 
results are expected to agree). 

The remarkable agreement between the results coming from both the BHM 
and the AS techniques contrast with those coming from the TSP probabilities. 



Therefore, this last one should be discarded as an approach to the simulation 
of physical systems in the TS framework. In Fig. 1 the values for q were chosen 
to be close to one. 

We are going to discuss separately the g < 1 and g > 1 regimes, since the 
system behaves differently in these regions. 



5.1 The q <1 regime 



Fig. 2 shows the internal energy as a function of T, for some values of q very 
close to one. The first surprising result is that a reentrant region develops as 
soon as q gets smaller than 1. This reentrant behaviour is even more noticeable 
as the lattice size increases, for the same value of q. This behaviour is already 
being reported [26] for a two-level system. We are going to use the same recipe 
proposed therein to deal with this pathology. 

The reentrant behaviour is also present in the free energy curve, as we can 
see in Fig. 3. To reconstruct the correct curve for the free energy and, as 
a consequence, for any thermodynamical macroscopic quantity, we have to 
choose a criterion to discard two of the three possible states in the loop (see 
the inset of Fig. 3). We choose to consider only the state that corresponds to 
the lowest value of the free energy. The correct curve for the internal energy is 
also displayed in Fig.3, and is another illustration that supports the statement 
that TS with renormalized q-expectation values is thermodinamically stable. 
A deeper discussion of the reentrant behaviour and the technique developed 
to restore uniqueness is being published elsewhere 





Fig. 2. Internal energy for L = 30 and L = 70. The curves (from left to right in 
all graphs) correspond to different q values ranging from 0.9990 to 1. The reentrant 
behaviour is being reported to appear in a two-level system. 

When uniqueness is restored by removing the loop, there results a free energy 
with a discontinuous first derivative with respect to temperature at a point 
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Fig. 3. Internal energy as a function of the temperature for L = 50 and q = 0.9990. 
The dashed hne shows the correct behaviour when non-uniqueness is removed by 
stabihty arguments. The inset shows the loop in the free energy that corresponds to 
the reentrant region, the removal of which restores uniqueness; the resulting internal 
energy is left with a discontinuous derivative with respect to temperature. 

which could be identified as the transition temperature. This means that the 
entropy has a discontinuity at this temperature. In addition, the magnetiza- 
tion, as shown in Fig. 4, also presents a discontinuity at this point - after the 
reentrancies have been removed. These discontinuities become more notice- 
able as we consider ever smaller values of q (see Fig. 4) and/or bigger lattices 
(see Fig. 5). All these results were obtained for finite lattices. To determine 
the transition temperature in the thermodynamic limit L — i> oo, we plot the 
transition temperature for finite lattices as a function of its inverse size, 1/L, 
and take the limit 1/L -^ 0. 
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Fig. 4. Curves obtained after the correction procedure discussed in the text. The 
plot of the internal energy on the left correspond to the right part of Fig. 2. On the 
right, we show the magnetization. The values of q are the same as in Fig. 2. 

Fig. 6 shows an attempt of finite-size scaling for different values of q. For q < 1 
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Fig. 5. Internal energy and magnetization as functions of temperature, now for a 
fixed value of q and different lattice sizes. 
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Fig. 6. Determination of the transition temperature. The lines are only guides to 
the eyes. 

and L ^ oo, the critical temperature vanishes extremely fast. Our results 
point out that there is no phase transition for g < 1 in the thermodynamic 
limit. 



5.2 The q > 1 regime 



Basically, the same approach we have used for g < 1 is also used here. However, 
we have found that the 2D Ising model presents a quite different behaviour in 
both regimes (actually, it is also different from the q = 1 case). Fig. 7 shows 
the internal energy and magnetization as functions of T for some values of q 
on a L = 70 lattice. From the internal energy curve, we promptly find that the 
discontinuity in its derivative (specific heat) at T = Tc is not present anymore. 
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therefore there are no evidence of a phase transition. This is also clear from the 
behaviour of the magnetization as a function of temperature. Fig. 8 supports 
this conclusion; the magnetization gets smoother with increasing values of q 
or the lattice size. 
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Fig. 7. Internal energy and magnetization as functions of temperature for L = 70 
and g = 1 to 1.0006. The boxes are used to display regions where abrupt changes 
in the derivative of the functions occur. The curves on this and the next figures are 
vertically displaced to allow better examination of these regions. 
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Fig. 8. Internal energy and magnetization as functions of temperature for a fixed 
value oi q = 1.0006 and for different lattice sizes. Again, the boxes are used to 
display regions where abrupt changes in the derivative of the functions occur. 

A careful examination of Fig. 7 and Fig. 8 reveals discontinuous changes in 
the derivative of the curves at two points (these changes are more easily seen 
for the largest values of q and V). These discontinuous derivatives are related 
to the transformation T -^ T' . Fig. 9 presents the relation between T and T' 
for the same values of L used in Fig. 8 and for a range of values of q slightly 
larger than what was used in Fig. 7. The curves display a region of abrupt 
change (but not a discontinuity) which gets larger as q is increased, for a 
fixed L. The values of T that limit this region are also those that lead to a 
discontinuous specific heat, and appear to have zero and oo as limits when q or 
L increase. In the thermodynamic limit, the magnetization is always positive 
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(therefore, the system is always in the ordered state, except at the hmit of 
infinite temperature) and the derivative of the specific heat is always positive. 
Therefore, we argue that there is no phase transition for the 2D Ising model, 
for q > 1. 



H 10 





Fig. 9. T as a function of T', for different lattice sizes. The curves correspond to 
different values of q ranging from 1 to 1.001. 



6 Conclusion 



We studied the simulation of magnetic systems in the Tsallis Statistics using 
the Broad Histogram Method. It was shown that this method is very effi- 
cient, since all thermodynamic observables of interest can be calculated for a 
new value of the parameter q without the need for a new computer run. The 
square lattice Ising model with nearest neighbour interactions was chosen as 
an example to test the method. All previous work on the nearest neighbour 
Ising model suggest that there is a phase transition for each value of g 7^ 1. 
However, our results show that there is no transition in this model for any 
value of g 7^ 1. Only the g = 1 system presents the usual second order phase 
transition. 

Further step along these lines would be the simulation of spin models with 
long-range interactions. We believe that these intrinsically non-extensive sys- 
tems can display a richer behaviour, within the Tsallis Statistics framework, 
than the model studied in this work. Actually, previous studies of magnetic 
systems in the Tsallis Statistics have already suggest that this is the case. 
However, a powerful simulational tool such as the Broad Histogram Monte 
Carlo Method used in this work was lacking and only recently has become 
available. The demonstration that we now present of its usefulness will surely 
bring considerable advance in the understanding of the behaviour of magnetic 
systems in a non-extensive regime. 
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